Target ecosystem

Target ecosystem is Feldmark, we use the shapefile with the final map:

arch <- sprintf("%s/all states feldmark/final outputs/all_states_feldmark_min.shp",MAPS)
eco.xy <- st_read(arch)
## Reading layer `all_states_feldmark_min' from data source `/srv/scratch/z3529065/gisdata/aust-alps/all states feldmark/final outputs/all_states_feldmark_min.shp' using driver `ESRI Shapefile'
## Simple feature collection with 237 features and 55 fields
## geometry type:  MULTIPOLYGON
## dimension:      XYZ
## bbox:           xmin: 384346.9 ymin: 5182875 xmax: 619083.8 ymax: 5970740
## z_range:        zmin: 0 zmax: 0
## CRS:            28355

Climate data is in coarse cells (10 km²), while ecosystems are very small (<40 ha)

eco.area <- st_area(eco.xy)
units(eco.area) <- with(ud_units, ha)
##units(eco.area) <- with(ud_units, km^2)
hist(eco.area)

Thus we can simplify the analysis using the centroids of the ecosystem polygons. We could use the areas as weights for the final summary of results if needed.

eco.centroids <- st_coordinates(st_transform(st_centroid(eco.xy),crs = "EPSG:4326"))
## Warning in st_centroid.sf(eco.xy): st_centroid assumes attributes are constant
## over geometries of x

Interactive maps with plotly:

map1 <- plot_mapbox(data.frame(eco.centroids)) %>% add_markers(x=~X,y=~Y,color = I("maroon"),size=10)

map1 %>% layout(
   mapbox = list(
      zoom = 8,
      center = list(lat = -36.5, lon = 148.5)
   )
)
## Warning: `arrange_()` is deprecated as of dplyr 0.7.0.
## Please use `arrange()` instead.
## See vignette('programming') for more help
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_warnings()` to see where this warning was generated.

Load GDD data for one model and scenario and aggregate GDD values per cell using inverse distance weighted interpolation:

if (!exists("gdd.eco")) {
   gdd.eco <- data.frame(eco.centroids)

   for (PERIOD in c("1990-2009","2020-2039","2060-2079")) {
      load(sprintf("%s/%s-%s-%s.rda",RDATA,PERIOD,MODEL,PRM))

      for (yy in unique(GDD$year)) {
         ss <- subset(GDD, year %in% yy & n>360 & lon>(min(eco.centroids[,1])-1) & lon<(max(eco.centroids[,1])+1) & lat>(min(eco.centroids[,2])-1) & lat<(max(eco.centroids[,2])+1))
         if(nrow(ss)>0) {

         dst1 <- pointDistance(ss[,2:1],eco.centroids,lonlat=T,allpairs=T)

         #Inverse distance weighting
         dst1[dst1>50000] <- NA
         w <- 1/dst1^3
         gdd.eco[yy] <- apply(w,2,function(x) sum(x*ss$GDD,na.rm=T)/sum(x,na.rm=T))

         }
      }   
   }   
}

valid.eco <- apply(!is.na(gdd.eco[,-(1:2)]),1,sum)>10
gdd.eco <- gdd.eco[valid.eco,]
head(gdd.eco)
##          X         Y     1990     1991     1992     1993     1994     1995
## 1 148.3284 -36.44514 1286.042 1283.316 1299.349 1361.629 1306.813 1387.313
## 2 148.3270 -36.44428 1285.492 1282.626 1299.068 1360.960 1305.906 1386.195
## 3 148.3287 -36.44484 1286.044 1283.341 1299.285 1361.645 1306.889 1387.443
## 4 148.2768 -36.43181 1256.750 1249.418 1272.203 1328.572 1272.414 1349.492
## 5 148.3038 -36.40235 1318.676 1319.917 1336.123 1395.673 1337.489 1416.236
## 6 148.2870 -36.41420 1304.245 1302.100 1321.250 1379.675 1321.959 1400.011
##       1996     1997     1998      1999     2000     2001     2002     2003
## 1 1455.817 1366.046 1281.125 1046.0798 1371.507 1307.054 1372.039 1175.881
## 2 1454.176 1365.946 1279.835 1044.2151 1370.128 1306.318 1371.423 1175.259
## 3 1456.035 1365.950 1281.277 1046.3707 1371.667 1307.084 1372.049 1175.901
## 4 1410.221 1339.921 1239.619  997.7693 1329.931 1274.149 1339.951 1144.194
## 5 1482.694 1406.644 1313.772 1071.6094 1397.744 1338.569 1407.711 1209.425
## 6 1464.625 1390.985 1294.162 1051.9924 1382.269 1323.381 1391.063 1193.453
##       2004     2005     2006     2007     2008     2020     2021     2022
## 1 1336.245 1443.444 1262.834 1342.791 1324.918 1553.624 1248.681 1727.750
## 2 1335.607 1442.570 1261.756 1341.959 1323.794 1552.522 1246.676 1727.547
## 3 1336.262 1443.509 1262.949 1342.846 1325.048 1553.752 1249.007 1727.660
## 4 1305.021 1407.202 1226.389 1307.872 1286.617 1515.709 1199.441 1700.493
## 5 1367.603 1479.055 1292.450 1376.140 1357.989 1586.722 1272.543 1767.706
## 6 1353.774 1461.310 1275.702 1359.525 1339.338 1567.834 1252.614 1752.661
##       2023     2024     2025     2026     2027     2028     2029     2030
## 1 1269.247 1366.308 1431.073 1590.592 1324.293 1581.255 1383.444 1555.677
## 2 1267.310 1365.320 1430.091 1590.201 1322.699 1580.625 1381.433 1555.562
## 3 1269.560 1366.412 1431.136 1590.537 1324.512 1581.260 1383.761 1555.562
## 4 1220.574 1328.953 1394.566 1560.517 1278.328 1547.505 1331.250 1529.362
## 5 1295.202 1403.032 1458.358 1628.528 1353.258 1621.025 1413.726 1595.772
## 6 1274.678 1383.182 1445.130 1612.996 1333.947 1602.214 1390.976 1580.248
##       2031     2032     2033     2034     2035     2036     2037     2038
## 1 1517.808 1223.788 1603.375 1521.348 1491.134 1523.302 1535.947 1636.687
## 2 1516.832 1221.711 1601.478 1519.127 1490.073 1522.401 1533.972 1636.079
## 3 1517.901 1224.128 1603.650 1521.689 1491.232 1523.347 1536.252 1636.677
## 4 1480.215 1171.862 1554.358 1463.959 1452.255 1486.240 1483.779 1602.936
## 5 1554.559 1251.822 1625.664 1552.825 1523.165 1555.894 1567.773 1674.423
## 6 1535.272 1229.811 1609.013 1528.412 1506.058 1540.024 1544.214 1657.585
##       2060     2061     2062     2063     2064     2065     2066     2067
## 1 1676.089 2090.012 1907.188 2061.659 1703.823 1857.639 1781.925 1996.550
## 2 1674.316 2088.733 1906.150 2059.869 1702.225 1857.048 1780.439 1995.503
## 3 1676.344 2090.142 1907.246 2061.902 1704.027 1857.634 1782.101 1996.643
## 4 1624.795 2046.657 1867.356 2013.509 1655.088 1823.869 1736.190 1954.565
## 5 1711.913 2122.775 1938.608 2086.893 1736.730 1897.323 1812.539 2038.497
## 6 1688.567 2105.230 1924.043 2070.079 1715.830 1879.347 1794.057 2016.347
##       2068     2069     2070     2071     2072     2073     2074     2075
## 1 1853.569 1960.387 1994.242 2126.437 1991.575 1833.775 1943.632 1898.518
## 2 1852.832 1959.306 1993.147 2125.682 1990.385 1833.466 1943.720 1897.719
## 3 1853.578 1960.477 1994.325 2126.435 1991.679 1833.688 1943.462 1898.535
## 4 1816.800 1918.656 1952.926 2087.845 1947.374 1803.348 1917.309 1859.001
## 5 1889.264 1999.857 2026.756 2166.440 2030.015 1872.355 1988.450 1941.340
## 6 1873.569 1979.155 2010.741 2147.933 2009.077 1857.629 1972.114 1920.868
##       2076     2077     2078
## 1 1921.658 2095.358 1904.635
## 2 1919.597 2094.340 1904.530
## 3 1921.947 2095.438 1904.509
## 4 1868.466 2055.480 1873.172
## 5 1946.205 2130.600 1954.237
## 6 1928.034 2113.427 1934.764

For one location, we can calculate the trend in annual GDD for the whole time period, and interpolate the expected values for years 2000 and 2050:

y <- unlist(gdd.eco[1,-(1:2)])
x <- as.numeric(colnames(gdd.eco)[-(1:2)])
mdl <- lm(y~x,subset=y>0)
IV <- predict(mdl,data.frame(x=2000))
FV <- predict(mdl,data.frame(x=2050))

plot(x,y,type='n',xlab="Year",ylab="GDD")
rect(2000,-3000,2050,3000,col="palegoldenrod",border="palegoldenrod")
points(x,y,col="maroon",pch=1.2)
abline(mdl,lty=2)
points(2000,IV,pch=19,cex=1.6,type="p")
points(2050,FV,pch=19,cex=1.6,type="p")
text(2000,IV*.9,sprintf("Initial value = %0.2f",IV),adj=c(0,1))
text(2050,FV*1.1,sprintf("Final value = %0.2f",FV),adj=c(1,1))

If we set an arbitrary collapse value of \(GDD[collapse]=2000\), the relative severity for this location is:

CT <- 2000

(FV-IV)/(CT-IV)
##         1 
## 0.6139087

We can now repeat this for all locations, and calculate the relative severity for all units:

x <- as.numeric(colnames(gdd.eco)[-(1:2)])
CT <- 2000
eco.RS <- data.frame()
for (k in 1:nrow(gdd.eco)) {
   y <- unlist(gdd.eco[k,-(1:2)])
   mdl <- lm(y~x,subset=y>0)
   IV <- predict(mdl,data.frame(x=2000))
   FV <- predict(mdl,data.frame(x=2050))
   eco.RS <- rbind(eco.RS,data.frame(k,RS=(FV-IV)/(CT-IV)))
}
summary(eco.RS$RS)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##  0.5715  0.5848  0.6169  0.6120  0.6294  0.6491

In this case, relative severity is > 50% for all localities (>80% extent), thus the resulting category would be EN.

We could use reference data to estimate a more appropriate collapse threshold. For example Hakea microcarpa was mentioned as a species that could invade Feldmark under rising temperatures. We use occurrence localities from the Atlas of Living Australia as of an external potentially invasive species)

arch <- sprintf("%s/Hakea-microcarpa-locs.csv",MAPS)
Hm.xy <- unique(read.table(arch,head=F))
coordinates(Hm.xy) <- 1:2
proj4string(Hm.xy) <-  '+proj=utm +zone=55 +south +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs '
Hm.ll <- spTransform(Hm.xy,"+init=epsg:4326")

The occurrence records are very widespread, in Tasmania the distribution of Hakea microcarpa has little overlap with the distribution of the ecosystem...

plot(Hm.ll,col=2,cex=.5)
points(eco.centroids,cex=.5)

But in the mainland the Feldmark has a very restricted distribution (2 cells from the NARCLiM domain), and it overlaps with records of H. microcarpa:

plot(Hm.ll,col=2,cex=.5,ylim=c(-37,-36),xlim=c(148,149))
points(unique(GDD[,2:1]),pch=22,cex=8)
points(eco.centroids,cex=.5)

We can calculate the expected GDD values for the locations with Hakea microcarpa:

if (!exists("gdd.col")) {
   xys2 <- coordinates(Hm.ll)
   gdd.col <- data.frame(xys2)

   for (PERIOD in c("1990-2009")) {
      load(sprintf("%s/%s-%s-%s.rda",RDATA,PERIOD,MODEL,PRM))

      for (yy in unique(GDD$year)) {
         ss <- subset(GDD, year %in% yy & n>360 & lon>(min(eco.centroids[,1])-1) & lon<(max(eco.centroids[,1])+1) & lat>(min(eco.centroids[,2])-1) & lat<(max(eco.centroids[,2])+1))
         if(nrow(ss)>0) {
            dst2 <- pointDistance(ss[,2:1],xys2,lonlat=T,allpairs=T)
            #Inverse distance weighting
            dst2[dst2>50000] <- NA
            w <- 1/dst2^3
            gdd.col[yy] <- apply(w,2,function(x) sum(x*ss$GDD,na.rm=T)/sum(x,na.rm=T))
         }
      }   
   }   
}
head(gdd.col)
##         V1        V2     1990     1991     1992     1993     1994     1995
## 1 149.7500 -36.20000 1702.095 1768.323 1735.654 1804.114 1795.474 1906.337
## 2 147.5648 -41.81434      NaN      NaN      NaN      NaN      NaN      NaN
## 3 148.7012 -35.73178 1210.436 1221.998 1236.690 1300.833 1267.680 1333.256
## 4 149.4308 -36.67100 1751.544 1828.651 1756.623 1832.247 1828.904 1970.016
## 5 146.7297 -41.98384      NaN      NaN      NaN      NaN      NaN      NaN
## 6 151.4442 -30.05397      NaN      NaN      NaN      NaN      NaN      NaN
##       1996     1997     1998     1999     2000     2001     2002     2003
## 1 2052.427 1817.797 1816.942 1561.766 1881.203 1837.704 1861.606 1685.009
## 2      NaN      NaN      NaN      NaN      NaN      NaN      NaN      NaN
## 3 1428.697 1344.327 1213.206 1027.469 1330.314 1269.503 1299.355 1147.693
## 4 2108.241 1840.853 1864.836 1610.293 1921.716 1863.949 1908.139 1715.518
## 5      NaN      NaN      NaN      NaN      NaN      NaN      NaN      NaN
## 6      NaN      NaN      NaN      NaN      NaN      NaN      NaN      NaN
##       2004     2005     2006     2007     2008
## 1 1787.898 1941.329 1734.138 1858.485 1787.209
## 2      NaN      NaN      NaN      NaN      NaN
## 3 1276.378 1396.949 1218.293 1289.275 1291.683
## 4 1812.001 1985.369 1778.215 1876.343 1829.545
## 5      NaN      NaN      NaN      NaN      NaN
## 6      NaN      NaN      NaN      NaN      NaN

But there is considerable overlap and no clear threshold to separate the ecosystem from the areas with presence of H. microcarpa

x <- as.numeric(colnames(gdd.col)[-(1:2)])
y <- as.numeric(colnames(gdd.eco)[-(1:2)])

boxplot(gdd.col[,-(1:2)],at=x,col=grey(.6),border=grey(.4))
matpoints(y,t(gdd.eco[,-(1:2)]),pch=1,col="darkgreen",cex=1)